//+------------------------------------------------------------------+
//|																			Mem2D.h	|
//| Lightweight (fully inline) templated matrix storage class.			|
//| Implemented as a one dimensional column major array,					|
//| which is, for example, LAPACK and OpenGL compatible.					|
//| I.e., incrementing the first index (row) will address				|
//| consecutive memory positions.												|
//| The beginning of the first element (0,0) is 64 byte aligned.		|
//| As a tradeof, data types must be 'simple', i.e.,						|
//| sizeof(T) must return the total memory of a T object.				|
//| A matrix with # columns n = 1, is an implicit column vector.		|
//+------------------------------------------------------------------+
#pragma once

//+------------------------------------------------------------------+
template<class T>
class Mem2D
{
private:
//////////////////////////////////////////////////////////////////////
	char*	data_unaligned;
	T*		data;
	int	rows;
	int	cols;

public:
//////////////////////////////////////////////////////////////////////
//+------------------------------------------------------------------+
	Mem2D()
	:	data_unaligned	(NULL),
		data				(NULL),
		rows				(0),
		cols				(0)
	{
	}

//+------------------------------------------------------------------+
	Mem2D(int m, int n = 1)
	:	data_unaligned	(NULL),
		data				(NULL),
		rows				(0),
		cols				(0)
	{
		resize(m,n);
	}

//+------------------------------------------------------------------+
	Mem2D(int m, int n, const T &init)
	:	data_unaligned	(NULL),
		data				(NULL),
		rows				(0),
		cols				(0)
	{
		resize(m,n);

		if(init)
			for(int i = 0; i < rows*cols; i++)
				data[i] = init;
		else
			zero();
	}

//+------------------------------------------------------------------+
	Mem2D(int m, int n, const T* init, bool transpose = false)
	:	data_unaligned	(NULL),
		data				(NULL),
		rows				(0),
		cols				(0)
	{
		resize(m,n);
		if(transpose)
			loadTranspose(init);
		else
			memcpy(data, init, sizeTotal());
	}

//+------------------------------------------------------------------+
//	copy contructor
	Mem2D(const Mem2D& M, bool transpose = false)
	:	data_unaligned	(NULL),
		data				(NULL),
		rows				(0),
		cols				(0)
	{
		if(transpose)
		{
			resize(M.cols,M.rows);
			for(int i = 0; i < rows; i++)
			for(int j = 0; j < cols; j++)
				(*this)(i,j) = M(j,i);
		}
		else
		{
		//	copy assignment operator
			*this = M;
		}

	}

//+------------------------------------------------------------------+
//	copy assignment operator
	Mem2D& operator = (const Mem2D& M)
	{
		if(this != &M)							// check for self-assignment
		{
			resize(M.rows,M.cols);
			memcpy(data, M.data, sizeTotal());
		}

		return *this;
	}

//+------------------------------------------------------------------+
//	destructor
	~Mem2D()
	{
		free();
	}

//+------------------------------------------------------------------+
//	The matrix will reference some other memory (e.g. a sub-matrix).
//	It does not (re)allocate or delete its own storage.
	Mem2D& reference(T* ref, int m = -1, int n = -1)
	{
		data = ref;
		if(m>=0)	rows = m;
		if(n>=0)	cols = n;
		return *this;
	}

//+------------------------------------------------------------------+
//	Will change the dimension of the Mem2D.
//	If possible, the memory is left untouched!
//	If copy is true and the size is increased or m is changed,
//	the old data is copied into the new memory, column by column.
	Mem2D& resize(int m, int n = 1, bool copy = false)
	{
		if(m<0)
			m=0;
		if(n<0)
			n=0;

	//----
	//	Need more memory?
		size_t bytes = m*n*sizeof(T);
		if(bytes > sizeReserved() || (copy&&rows!=m))
		{
		// Can we free the old memory _before_ allocating new?
			if(!copy)
				delete[] data_unaligned;

			char	*new_data_unaligned = new char[bytes+96];
			T		*old_data = data;

		//----
		//	Make a 64 byte aligned pointer
			data = (T*)(((uintptr_t)(new_data_unaligned+96) >> 6) << 6);

			if(copy)
			{
				if(rows!=m)
				{
				//	need to copy each column separately
					int r = min(rows, m);
					for(int j = 0; j < cols; j++)
						memcpy(data+j*m, old_data+j*rows,  r*sizeof(T));
				}
				else
						memcpy(data,     old_data, cols*rows*sizeof(T));

				delete[] data_unaligned;
			}

			data_unaligned = new_data_unaligned;

		// save size of (aligned) memory
			*(size_t*)data_unaligned = bytes;
		}

	//----
		rows = m;
		cols = n;
		return *this;
	}

//+------------------------------------------------------------------+
//	Does not touch the memory, so make sure that it's big enough!
	Mem2D& assumeSize(int m, int n = 1)
	{
		rows = m;
		cols = n;
		return *this;
	}

//+------------------------------------------------------------------+
	size_t sizeTotal() const
	{
		return rows*cols*sizeof(T);
	}

//+------------------------------------------------------------------+
	size_t sizeReserved() const
	{
		return data_unaligned ? *(size_t*)data_unaligned : 0;
	}

//+------------------------------------------------------------------+
//	Set every matrix element to zero
	Mem2D& zero()
	{
		memset(data,0,sizeTotal());
		return *this;
	}

//+------------------------------------------------------------------+
//	Delete the allocated memory and set matrix size to 0
	Mem2D& free()
	{
		delete[] data_unaligned;
		data_unaligned	= NULL;
		data				= NULL;
		rows				= 0;
		cols				= 0;
		return *this;
	}

//+------------------------------------------------------------------+
	Mem2D& initialize(const T &x)
	{
		T	*dst = data,
			*lim = data+rows*cols;

		while(dst < lim)
			*dst++ = x;

		return *this;
	}

//+------------------------------------------------------------------+
	int m() const	{	return rows;	}
	int n() const	{	return cols;	}

//+------------------------------------------------------------------+
//	Will increment the number of rows and copy one element.
	T& push_back(const T &x)
	{
		int i = rows;
		resize(rows+1,cols,true);
		return data[i] = x;
	}

//+------------------------------------------------------------------+
//	Will return the last element of the first column.
	const	T&	back()								const	{	return data[rows-1];		}
			T&	back()										{	return data[rows-1];		}
	const	T&	operator()	(int i, int j)	const	{	return data[i+j*rows];	}
			T&	operator()	(int i, int j)			{	return data[i+j*rows];	}
	const	T&	get			(int i, int j)	const	{	return data[i+j*rows];	}
			T&	get			(int i, int j)			{	return data[i+j*rows];	}
	const	T*	getCol		(int j)				const	{	return data + j*rows;	}
			T*	getCol		(int j)						{	return data + j*rows;	}

//+------------------------------------------------------------------+
//	conversion operator returns the single dimensional data array,
//	i.e. the address of the first element.
	operator	const	T*		()	const						{ return data; }
	operator			T*		()								{ return data; }
	operator	const	void*	()	const						{ return data; }
	operator			void*	()								{ return data; }

//+------------------------------------------------------------------+
	Mem2D<T> &transpose()
	{
		if(rows==cols)							//	In place transformation for square matrices
		{
			for(int i = 0;   i < rows; i++)
			for(int j = i+1; j < cols; j++)
				swp((*this)(i,j), (*this)(j,i));
		}
		else										//	Temporary copy for rectangular matrices (~5-10x slower)
		{
			const Mem2D tmpTranspose(*this,true);
			*this = tmpTranspose;
		}

		return *this;
	}

//+------------------------------------------------------------------+
	Mem2D<T> &loadTranspose(const T* init)
	{
		for(int i = 0; i < rows; i++)
		for(int j = 0; j < cols; j++)
			data[i+j*rows] = init[i*cols+j];

		return *this;
	}

//+------------------------------------------------------------------+
// only load one source column
	Mem2D<T> &loadTranspose(const T* init, int col)
	{
		for(int j = 0; j < cols; j++)
			data[col+j*rows] = init[col*cols+j];

		return *this;
	}

//+------------------------------------------------------------------+
//	Will compare each element to x and return a pointer to the element,
//	if found, or NULL if not.
	T *find(const T &x)
	{
		for(int i = 0; i < rows*cols; i++)
			if(data[i] == x)
				return &data[i];

		return NULL;
	}

};//class Mem2D



#include <ostream>
//+------------------------------------------------------------------+
//	human readable output
template<class T>
inline std::ostream& operator << (std::ostream &out, const Mem2D<T> &M)
{
	for(int i = 0; i < M.m(); i++)
	{
		out << "[";
		for(int j = 0; j < M.n(); j++)
			out << M(i,j) << (j<M.n()-1 ? " " : "]\n");
	}

	return out;
}




//+------------------------------------------------------------------+
//	Has an additional index, to keep track of the next position.
template<class T>
class CircleBuffer : public Mem2D<T>
{
	int index;

public:
//+------------------------------------------------------------------+
	CircleBuffer(int size = 0)
	:	Mem2D<T>(size),
		index		(0)
	{
		zero();
	}

//+------------------------------------------------------------------+
//	Will copy the element to the position of the pointer,
//	increment the pointer and reset to 0, if out of bounds.
	CircleBuffer& save(const T &x)
	{
		operator T*()[index] = x;
		if(++index == m())
			index = 0;
		return *this;
	}

};//class CircleBuffer


